///\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
//(c) Tashan and Friends                 \\
//Generalised 1D Riemann Problem Solver  //
//                                       \\
//\\//\\//\\//\\/\\//\\//\\//\\//\\/\\/\\//


//Just change the parameters in init functions, dont chance the size of the domain, however even that can be changed with slight modification.

#include <iostream>
#include <fstream>
#include <time.h>
#include <sys/time.h>
#include <stdlib.h>
#include <cmath>
#include <string>
#define pi 3.1415926536
using namespace std;



class Reimann{
public:
  //Variables Description
  //
  //size: size of grid
  //t: time at which solution to be calculated
  //Ur,Ul: right and left values of properties at t=0, U={rho,p,u,T,a}
  //gamma ratio of specific heat, R ideal gas constant
  //Ar,Al,Br,Bl coeff in pressure equation
  //pstar: solution of the f(p), pinit: initial value used in solving f(p)=0
  //ustar: value of u where p=pstar
  int size;
  float t;
  double rhol,rhor,pl,pr,ul,ur,Tl,Tr,al,ar;
  double Al,Ar,Bl,Br;
  float gamma, R;
  double pstar, pinit;
  double ustar; 
  char _filename[100];
  
  //Functions
  //_init: Assign values to the variables
  //_calcAB: Calculate AL and BL
  //_getfunc(p): return the value of function of f(p)
  //_getder(p): return derivative of the function f(p)
  //_initp: calculate initial value of p to be used in calculating root
  //_getroot(p): gives root of function f(p) using newton raphson
  //_getu: calculate ustar based on p
  //_getdata: sample the points and assingned the appropriate value of properties
  //_writedata: write data to file
  //letsgo: calls all the function in sequence to slolve the problem for 
  //        the given initial condition
  void _init(void);
  void _calcAB(void);
  double _getfunc(double p);
  double _getder(double p);
  void _initp(void);
  double _getroot(double p);
  void _getu(void);
  void _getdata(void);
  void _writedata(void);
  void letsgo(void);
  //Geometry
  double x[100];
  double p[100];
  double rho[100];
  double u[100];
  double E[100];

};
  
void Reimann::_init(void){
  sprintf(_filename,"test_4.txt");
  gamma = 1.4;
  t = .035;

  //Left State
  rhol = 1.0;
  ul = 0.0;
  pl = 0.01;
  //Right State
  rhor = 1.0;
  ur = 0.0;
  pr = 100.0;
    
  al = sqrt(gamma*pl/rhol);
  ar = sqrt(gamma*pr/rhor);
  
  for(int i=0;i<100;i++){
    x[i] = -0.5+i*.01;
  }

}

void Reimann::_calcAB(void){
  Al = 2/((gamma+1)*rhol);
  Ar = 2/((gamma+1)*rhor);
  Bl = (gamma-1)*pl/(gamma+1);
  Br = (gamma-1)*pr/(gamma+1); 
}

double Reimann::_getfunc(double p){
  double func = 0.0;
  if (p>pl){
    //Shock Wave
    func += (p-pl)*sqrt(Al/(p+Bl));
  }
  else{
    //Rarefraction Wave
    func += (2*al/(gamma-1))*(pow(p/pl,(gamma-1)/2/gamma)-1);
  }
  if (p>pr){
    //Shock Wave
    func += (p-pr)*sqrt(Ar/(p+Br));
  }
  else{
    //Rarefraction Wave
    func += (2*ar/(gamma-1))*(pow(p/pr,(gamma-1)/2/gamma)-1);
  }
  //Delta u
  func += (ur-ul);
  return func;
}
double Reimann::_getder(double p){
  double func = 0.0;
  if (p>pl){
    func += sqrt(Al/(Bl+p))*(1- (p-pl)/(2*(Bl+p)));
  }
  else{
    func += (1/(rhol*al))*pow(p/pl,-(gamma+1)/(2*gamma));
    }
  if (p>pr){
    func += sqrt(Ar/(Br+p))*(1- (p-pr)/(2*(Br+p)));
  }
  else{
    func += (1/(rhor*ar))*pow(p/pr,-(gamma+1)/(2*gamma));
    }
  
  return func;
}




double Reimann::_getroot(double p){
  double tmp;
  double error = _getfunc(p);
  while(abs(error)>1e-6){
    tmp = _getfunc(p)/_getder(p);
    p = p - tmp;
    error = _getfunc(p);
  }
  return p;
}

void Reimann::_initp(void){
  pinit = al+ar-.5*(gamma-1)*(ur-ul);
  pinit /= al/pow(pl,(gamma-1)/2/gamma) + ar/pow(pr,(gamma-1)/2/gamma);
  pinit = pow(pinit,2*gamma/(gamma-1));
  
}

void Reimann::_getu(void){
  double fr,fl;
  if (pstar>pl){
    //Shock Wave
    fl = (pstar-pl)*sqrt(Al/(pstar+Bl));
  }
  else{
    //Rarefraction Wave
    fl = (2*al/(gamma-1))*(pow(pstar/pl,(gamma-1)/2/gamma)-1);
  }
  if (pstar>pr){
    //Shock Wave
    fr = (pstar-pr)*sqrt(Ar/(pstar+Br));
  }
  else{
    //Rarefraction Wave
    fr = (2*ar/(gamma-1))*(pow(pstar/pr,(gamma-1)/2/gamma)-1);
  }

  ustar = .5*(ul+ur+fr-fl);

}

void Reimann::_getdata(void){
  double g = (gamma-1)/(gamma+1);
  double rhoo,S,Sh,St,astar;
  double rat;
  for(int i=0;i<100;i++){
    rat = x[i]/t;
    
    if(rat <= ustar){
     
      
      //Left of contact wave
      if(pstar>pl){
	//Shock
	rhoo = rhol*(pstar/pl+g)/(g*pstar/pl+1);
	S = ul - al*sqrt((gamma+1)*pstar/pl/2/gamma + (gamma-1)/2/gamma);
	if(x[i]/t<=S){
	  //Left State
	  rho[i] = rhol;
	  p[i] = pl;
	  u[i] = ul;
	  
	}
	else{
	  //Left Star state
	  rho[i] = rhoo;
	  p[i] = pstar;
	  u[i] = ustar;
	}
  
      }
      else{
	//Rarefraction
	rhoo = rhol*pow(pstar/pl,1/gamma);
	astar = al*pow(pstar/pl,(gamma-1)/2/gamma);
	Sh = ul-al;
	St = ustar - astar;

	if(rat<=Sh){
	  //Left State
	  rho[i] = rhol;
	  p[i] = pl;
	  u[i] = ul;
	
	}
	else if(rat>=St){
	  //Left Star State
	  rho[i] = rhoo;
	  p[i] = pstar;
	  u[i] = ustar;
	
	}
	else{
	
	  //Inside rarefraction wave
	  rho[i] = rhol*pow(2/(gamma+1) + g/al*(ul-rat),2/(gamma-1));
	  u[i] = 2/(gamma+1)*(al + (gamma-1)/2*ul + rat);
	  p[i] =pl*pow(2/(gamma+1) +g/al*(ul-rat),2*gamma/(gamma-1));
	    

	}
      }
    }
    else{
     
      //Right of contact wave
      if(pstar>pr){
	//Shock
	rhoo = rhor*(pstar/pr+g)/(g*pstar/pr+1);
	S = ur + ar*sqrt((gamma+1)*pstar/pr/2/gamma + (gamma-1)/2/gamma);
	if(rat>=S){
	  //Right State
	  rho[i] = rhor;
	  p[i] = pr;
	  u[i] = ur;
	  
	}
	else{
	  //Right Star state
	  rho[i] = rhoo;
	  p[i] = pstar;
	  u[i] = ustar;
	}
      }
      else{
	//Rarefraction
	rhoo = rhor*pow(pstar/pr,1/gamma);
	astar = ar*pow(pstar/pr,(gamma-1)/2/gamma);
	Sh = ur+ar;
	St = ustar + astar;
	  
	if(rat>=Sh){
	  //Right State
	  rho[i] = rhor;
	  p[i] = pr;
	  u[i] = ur;
	}
	else if(rat<=St){
	  //Right Star State
	  rho[i] = rhoo;
	  p[i] = pstar;
	  u[i] = ustar;
	}
	else{
	  //Inside rarefraction wave
	  rho[i] = rhor*pow(2/(gamma+1) - g/ar*(ur-rat),2/(gamma-1));
	  u[i] = 2/(gamma+1)*(-ar + (gamma-1)/2*ur + rat);
	  p[i] =pr*pow(2/(gamma+1) - (gamma-1)/(gamma+1)/ar*(ur-rat),2*gamma/(gamma-1));
	    

	}
	
      }
      
    }
    E[i] = p[i]/(gamma-1)/rho[i];
  }
}

void Reimann::_writedata(void){
   ofstream file (_filename);
   for(int i=0;i<100;i++){
     file<<x[i]<<" "<<p[i]<<" "<<rho[i]<<" "<<u[i]<<" "<<E[i]<<"\n";
   }
}

void Reimann::letsgo(void){
  _init();
  _calcAB();
  _initp();
  pstar = _getroot(pinit);
  _getu();
  _getdata();
  _writedata();
}

int main(){
  Reimann a;
  a.letsgo();
  return 0;

}
